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Abstract 

One of the methods for studying the highest energy cosmic rays is to measure the fluorescence light emitted by the 
extensive air showers induced by them. To reconstruct a shower cascade curve from measurements of the number of 
photons arriving from the subsequent shower track elements it is necessary to take into account the multiple scatterings 
that photons undergo on their way from the shower to the detector. In contrast to the earlier Monte-Carlo work, we 
present here an analytical method to treat the Rayleigh and Mie scatterings in the atmosphere. The method consists 
in considering separately the consecutive 'generations' of the scattered light. Starting with a point light source in 
a uniform medium, we then examine a source in a real atmosphere and finally - a moving source (shower) in it. 
We calculate the angular distributions of the scattered light superimposed on the not scattered light registered from 
a shower at a given time. The analytical solutions (although approximate) show how the exact numerical results 
should be parametrised what we do for the first two generations (the contribution of the higher ones being small). 
Not allowing for the considered effect may lead to an overestimation of shower primary energy by ~ 15% and to an 
underestimation of the primary particle mass. 
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1. Introduction 

One of the methods for studying extensive air show- 
ers of high energies (> 10 17 eV) is to register their im- 
ages in the optical (mainly fluorescence) light. This can 
be done by observing showers from the side in order to 
avoid the more intense Cherenkov light emitted rougly 
in the shower direction. The observations are made by 
a number of optical telescopes, each containing a large 
mirror and a camera with a matrix of photomultipliers 
(PMTs) placed at the focus of the optical system (HiRes 
yj], The Pierre Auger Observatory [2], The Telescope 
Array so that photons arriving from a given direc- 
tion on the sky are focused on a particular PMT (pixel). 
Photon arrival time can also be measured if the time 
structure of the PMT signals is recorded. A cosmic ray 
induced shower produces at a given time a light spot 
on the camera which moves across it as the shower de- 
velops in the atmosphere so that succeeding PMTs are 
being hit. 

It would have been ideal if the light producing a 
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shower image had contained the fluorescence photons 
only. This is because there exists experimental evidence 
that the number of fluorescence photons induced by a 
charged electron in the atmosphere is proportional to 
the energy lost by it for ionisation |4]. As practically 
all primary particle energy is eventually used for ioni- 
sation, this energy can be determined by measuring the 
fluorescence light emitted along the shower track in the 
atmosphere. 

However, there are several problems in deriving the flux 
of the fluorescence light emitted by a shower from that 
arriving at the detector. Firstly, the arriving light con- 
tains not only the fluorescence but also Cherenkov pho- 
tons. If the viewing angle (the angle between the line of 
sight and the shower direction) is large (say, > 30°) then 
it is mainly the Cherenkov light scattered in the atmo- 
sphere region just passed by the shower and observed 
by the detector (this light, before being scattered, trav- 
els roughly along the directions of the shower particles). 
Typically its fraction at the detector is about 15% of the 
fluorescence flux. For smaller viewing angles it is the 
Cherenkov light produced at the observed part of the 
shower that may dominate even the fluorescence signal. 

January 20, 2012 



The contribution of the Cherenkov light, which has to 
be subtracted from the total signal, has been extensively 
studied IS la III ISIl • 

The subject of this paper is another phenomenon, af- 
fecting shower images, most commonly called the mul- 
tiple scattering (MS) of light. Photons produced at the 
observed shower element, whatever their origin (fluo- 
rescence or Cherenkov), may undergo scattering in the 
atmosphere on their way from the shower to the detec- 
tor, causing an attenuation of the light flux arriving at 
the detector and a smearing of the image. This scat- 
tering may take place on the air molecules (Rayleigh 
scattering) or on larger transparent particles, aerosols 
(Mie scattering). Most of the scattered photons change 
their directions, so that they no longer arrive at the pixel 
registering the not scattered (direct) light. Moreover, 
they arrive later having longer path lengths to pass. On 
the other hand, photons emitted by the shower at ear- 
lier times and scattered somewhere, may fall in the field 
of view of the pixels just registering the direct photons 
emitted at a later time. The net effect is that the scattered 
light forms its own instantaneous image superimposed 
on that in the direct light. 

Our aim is to calculate the shower images in the multi- 
ply scattered light, so that this effect could be allowed 
for (subtracted) when determining the shower primary 
energy from the PMT signals. This problem was already 
studied by Roberts [9n, by us in several, short confer- 
ence contributions II 1 Oil and more recently by Pekala et 
al 111 111 . The approach of the other authors was based 
on Monte Carlo simulations of photons emitted by a 
shower. Photons were followed up to 5-6 scatterings 
and their arrival directions and time were registered by 
the detector. Many shower simulations were needed to 
obtain the MS images for various distances, heights, 
viewing angles of the observed shower parts. Finally, 
a phenomenological parametrisation of the number of 
MS photons was made as a function of the parameters 
found as relevant. 

In contrast to Roberts and Pekala et al this approach is 
based on an analytical treatment. The main idea is to 
consider the arriving MS light as a sum of the photons 
scattered only once (the first generation), of those scat- 
tered two times (the second generation) and so on and 
calculate separately the angular and temporal distribu- 
tions for each generation. 

We start (Section 2) with a consideration of the simplest 
situation when a point source of isotropic light flashes 
for a very short time in an uniform medium. We de- 
rive analytical expressions for the angular and temporal 
distributions of the first and next generations of light ar- 
riving at a particular distance from the source. 



As our aim is to apply our results to cosmic ray showers 
we need to consider a non-uniform medium like the at- 
mosphere. Assuming an exponential distribution of the 
gas density and similarly for aerosols we show that an 
effective scattering length between any two points in the 
atmosphere can be easily calculated analytically. Sig- 
nals of the first two generations arriving at a particular 
detector within a given angle £ to the direction to the 
source are found as a function of time (Section 3). 
Using these it is straightforward to derive the corre- 
sponding distributions if the source moves across the 
atmosphere, integrating the point source distributions 
over changing distance and time of light emission. In 
Section 4 we consider a moving light source, modelling 
a distant shower. We calculate angular distributions of 
MS light arriving at a detector at the same time as the di- 
rect (not scattered) photons emitted by the shower. This 
particular approach is quite natural because the data 
from optical detectors consist of the recorded signals by 
the camera PMTs within short time intervals At so that 
one needs to know how much of the MS light has to be 
subtracted from the main, direct signal. The method for 
calculating images in the MS light simultaneous with 
images in the direct light is relatively simple (for the 
first and the second generations) as it is based on the ge- 
ometry of the scattered photons in the particular gener- 
ation. It does not require time-consuming Monte-Carlo 
simulations that were done for various shower-observer 
geometries. Making some approximations we derived 
analytical formula for a shower image produced by the 
first generation (Section 4.2). Our analytical derivations 
allowed us to choose easily the variables on which and 
how to parametrise the MS signals. They also made us 
to realise that the dependence of it on the viewing angle 
was different for Rayleigh and Mie scatterings. Thus, 
we introduced a new, simple parametrisation of the frac- 
tion of the scattered photons arriving at the telescope 
within a given viewing cone, depending on the viewing 
angle of the shower (what has not been done before), 
separately for the two scatterings. For the same reason 
we also parametrised separately the second generation 
(Section 4.3). 

A discussion of the results and the implication of the MS 
effect on the derivation of shower parameters is given in 
Section 5. The last Section (6) contains a summary and 
conclusions. 

2. Point source flashing isotropically in uniform 
medium 

At any fixed time a distant shower can be regarded as 
a point source emitting isotropically fluorescence light 
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(about Cherenkov light see later). As explained above, 
we treat the light scattered in the medium as a sum 
of consecutive generations consisting of photons scat- 
tered only once, twice, and so on, on their way from the 
source to the observation point. 




Figure 1 : Geometry of the first scattering in a uniform medium. The 
light source is at the centre O of a sphere with radius R. Two rays 
shown are scattered at points Si and S2, correspondingly, arriving at 
the surface of the sphere with radius R at angles 9\ and 62. 
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where t = ct/R and c is the speed of light. 
The Jacobian of the transformation gives 
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Finally, the number of photons arriving at a unit surface 
at an angle (6, 8 + dff) (all azimuths) at time (f, t + dt) 
equals 
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2.1.1. Rayleigh scattering 

For the Rayleigh scattering we have 



2.1. First generation 

Let us consider the first generation, consisting of pho- 
tons scattered once only. We shall calculate the flux 
of these photons, j\{6,t;R), at a distance R from the 
source, such that j 1 (9,t;R)dCldtdS ± is the number of 
photons scattered only once, arriving at time (f, t + dt) 
after the flash, within a solid angle dQ.{0) at the surface 
dS j_ (perpendicular to the arrival direction) located at a 
distance R. To do this we shall calculate first the num- 
ber of photons crossing the sphere of radius R (from in- 
side) at an angle (6, 6+d6) with respect to the normal, at 
time (t, t + dt); see Fig. 1. The average number of pho- 
tons, per one photon emitted, interacting at a distance 
(x, x + dx) from the source and scattered at an angle a 
within d£l(a), equals 



_x dx 

dni(x, a) — e A — f(a)d£l(a) ■ e 
A 



(1) 



where A is the mean scattering path length, f(a)dQ. is 
the probability that, once the scattering has occured, 
the scattering angle is a within dQ.(a) = 2n sina da. To 
each pair of variables (x, a) there corresponds another 
pair (6, t) related to the former by 



a t — cosO 
t8 2 ~ sin6 



(2) 



f(a)=f R (a) = — (l+cos 2 a) 
16n 



(8) 



Expressing a as a function of and r (Eq. 2) we obtain 
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where y — r 2 - 2tcos0 + 1. Thus, the flux j^(9, t; R) of 
the first generation, defined above, equals 



3c e * 
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(10) 



where A R is the mean free path length for the Rayleigh 
scattering. However, A in the exponent depends on all 
the scattering processes active. In general it is deter- 
mined by 



1 = vi 



(11) 



for n processes. Thus, if both molecular and aerosol 
scatterings are active but one wants to calculate the flux 
of photons scattered by the Rayleigh process only, A in 
the exponent equals A-{ j- R + j^T 1 t> ut i n tne denomi- 
nator one has A R . 
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From Eq. 10 one can find the number of photons 

dNHr.CR) . . . . „ ... 

— ^ — arriving per unit time at a unit surface within 
a given angle as a function of time. 
We have 



dNf(t; £,R) 
dt 
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The integral can be found analytically, giving the result 
dNf(t; £,K) 



dt 
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where yi = (t — l) 2 , y% — t 2 - 2t cosf + 1, 02 = 1, 
fll = -4, fl = 6(t 4 + 1) - 4r 2 , a_i = -4(t 2 - l) 2 , 

fl_2 = (T 2 - l) 4 

We have also calculated analytically a similar distri- 
bution dN" /dt if the scattering was isotropic, i.e. if 
/(or) = i (Appendix A). 

One can also find analytically the angular distribution 

dM K 

of the arriving light (integrated over time), but for 
small angles only (Appendix B). The result is 
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where «: 1, k R = ~, if there is no Mie scattering and 
Ceu — 0.577 is the Euler constant. 
The ratio of all photons arriving within a small angle 
to those not scattered No, equals 

~ fl dMf . . 47rfl 2 9 k R e- k " 



1 r dM* 

— -dff 

No J dff 



e -k R 64 R 2 



(15) 
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where terms ~ 2 have been neglected. 

2.1.2. Mie scattering 

In the next paragraph we shall consider the scattering 
of light emitted by showers developing in the real at- 
mosphere, i.e. with the density depending on height. 



In addition to the Rayleigh process one has to take 
into account the Mie scattering occurring on particles 
(aerosols) larger that the light wavelength. The Mie an- 
gular distribution is concentrated at rather small angles, 
in contrast to the Rayleigh case. Moreover, in the deeper 
parts of the atmosphere the mean free path length for the 
Mie scattering may be comparable to that for Rayleigh, 
so that it is necessary to calculate the distribution of the 
light scattered by the Mie process only. 
As before, we start with a simpler case - a uniform 
medium. The angular distribution of light scattered on 
particles with sizes larger than the light wavelength de- 
pends on the distribution of the sizes and is not a well 
known function. Roberts [9] adopts a function of the 
form 



f(a) 



-Ba 



+ Ce 



Da 



(16) 



Here, however, we prefer an expression allowing us to 
perform some integrations analytically. Most crucial is 
to have the number of numerical integrations for the sec- 
ond generation as few as possible. We shall see that to 
find jj{0, t; R) for the Rayleigh scattering (Section 2.3) 
there is only one integration (over x') to be done numer- 
ically since the form of f R (a') enables one to integrate 
analytically over (/>' and 9' (Eq. 25 and 26). Thus, we 
adopt the following form for the Mie angular distribu- 
tion: 



ff(a) = a 2 cos"a + b for -<a<n (17) 

where a x = 0.857, a 2 = 0.125, b = 0.025. 
This function is normalised as follows 



a\cos a + b for 



IT 

< a < - 
2 

n 



f 

Jo 



f(a) ■ 2n sina da — 1 



(18) 



It describes quite reasonably the distribution used by 
Roberts. 

From Eq. 6 and 7 we have 

cr- kT 
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where i = 1 if tgj = - 
and i = 2 if tgj > 1 
Since 
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we obtain 
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(20) 
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where y — t 2 - 2t cos8+l, and ;' is determined as before . 
In principle, it is possible to find analytically the number 

„ , dNf(t;C,R) . . ... , „ , 

of photons — l -j t — arriving within an angle £ after time 
t per unit time. However, each of the nine integrals 



sin 6 

y 



sinOcosOdB 



(22) 



contains many terms itself, so that an analytical depen- 
dence on £ and/or t would be practically lost. Thus, we 

dN M 

have found -jj- by integrating the flux (Eq. 21) numer- 
ically (similarly to Eq. 12). 

Finally, the total flux of the first generation is the sum of 
the two fluxes arising from the two active mechanisms 
of the scattering 



2.2. The second generation 



(23) 




Figure 2: Geometry of the second and higher scatterings in a uniform 
medium. The light source is at O and the detector is the surface of the 
sphere with radius R. The picture shows the last scattering. 

These are the photons scattered exactly two times. 
We shall consider first a general case when there are 
more than one scattering processes (as Rayleigh and 
Mie). Let us call the process of the first scattering as 
A and this of the second one as B. Both A and B can be 
either Rayleigh or Mie. We denote their mean scattering 
path lengths by A A and A B , and the angular distribution 
functions of the scattering by f A (a) and f B a), corre- 
spondingly. The light source flashes isotropically at the 
centre of a sphere with radius R at time t = (Fig. 2). 
As before we want to calculate the number of photons 
crossing the surface of the sphere from inside at a given 
angle 6, at time f, per unit time. 

Let us consider the photons scattered for the second 
time at a distance R' from the source. The number of 
photons incident on a small surface AS ' at an angle 6' 
(within dD!) at time (f, t' + dt') and scattered within a 



spherical shell of thickness dR' by an angle a' (within 
d£l(a')) equals 



dx 

j A {0, <\ R')dQ'dt'AS'\cosd'\—f B (a')dCl 



(24) 



where dx = dR'/\cos6'\. As now both processes are 
active the meaning of A in the factor in the expres- 
sion for j A is the effective mean path length for both 
processes. 

The direction of the scattered photons is at an angle a 
to the radius of the sphere and c/Q = sina da d<p, where 
<p is the azimuth of the photons scattered for the second 
time. For any given direction (9',(f>') before and (a, 0) 
after the second scattering, the scattering angle a' fulfils 
the relation co sa' = cosacosQ' - sina sinff cos<p'. The 
only function depending on the azimuth angle (/>' of the 
incident photons is f B (a'). Denoting 



F B (9',a)= C f{a'W 
Jo 



(25) 



and integrating (25) over 6' we obtain for the number of 
photons incident on AS ' and scattered within dR' into 
the solid angle dQ.(a) the following expression 

dt'AS'dR'dQ. [ j A (8',t';R')F B (6',a)sin8'd6' = 
Jo 

= G AB (R', a) ■ dt'AS'dR'dQ. (26) 

where the function G AB (R', a) is defined by the integral 
in the l.h.s. of (26). 

The pair of fixed variables R' and a defines another pair 
x 1 and 6, where x' is the photon path length after the 
second scattering. The Jacobian of the transformation 
gives the relation 



R 

dR' sina da = —sin8\cos6\dx'd8 



(27) 



Putting AS ' = 4-nR' 2 , the contribution of photons scat- 
tered for the second time at a distance (R',R' + dR') to 
arrive at an angle (6, 8 + d8) at the sphere with radius R 
equals 



dn AB (0,t,x';R) = e-- G 



Mi 



R'(x',8),a(x',8) 



■4nR 2 ■ 2nsin8\cos8\ dx'dddt 



(28) 



The factor e~ x ' /A multiplied by e~ c, ' IA in the expression 
for j A gives e~ c 'l A , independent of x'. Integration over 
x 1 gives the total number of the above photons 
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Figure 3: Comparison of the first (dN\ldr) (three upper curves) and the second (dN2ldr) (three lower curves) generations as functions of time 
(e = T - 1 = ct/R - 1). Number of photons are within angle (. Uniform medium, R = 1. a). Rayleigh (solid lines) and isotropic (dashed lines) 
scattering, b). Two scattering processes at work: Rayleigh and Mie, each with A = 2R. 



where G AB = G AB /e- cl ' /A . 

The maximum value of x' results from fixing time t . We 
have that 

r 2 -l 



R 



x max ~ ct R ~ 2 



T - COSO 

Thus, the flux of the second generation AB equals 

•AS _ 1 

h 



d 2 n AB 



4nR 2 ■ 2nsin6\cos9\ d9dt 



(30) 



(31) 
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Gt B dx' 



and, integrated over 6 for 9 < gives dN AB (t; £,R)/dt. 
With the Rayleigh and Mie processes active we must 
take into account all four cases A = R or M and B — R or 
M. Finally, the flux of the combined second generation 
photons is a sum of all specific fluxes 



;RR . :RM . ;MR . MM 

J2 = J 2 + h + h + h 



(32) 



Some of the integrals defined in this Section can be 
found as analytical functions (Appendix C). It is of 
some importance when calculating higher generations 
(see the next Section). 



2.3. The next generations 

Any next generation of the scattered photons can be 
calculated in the same way as the second one has been 
found from the previous one (the first). To calculate the 
flux ji(6, t; R) of the i - th generation, given _/',■_ i (6, t; R) 
we proceed as before when calculating the second gen- 
eration from the first one (Eq. 24). The number of 
photons, incident on AS' at an angle 6',<p' within d£l' 
at time (t',t' + dt') and scattered along dx into dCl(a) 
equals: 



dx 

ji- 1 (0 , t' ■ R')dt'dQ.'AS '\cos0 1 —fi(a')dQ 
If there are two scattering processes, R + M, then 
Ma') 



A, 



f R {a') | f M {a') 



(33) 



(34) 



is the scattering probability by an angle a' by any pro- 
cess per unit distance per unit solid angle. The rest of 
the derivation of jj(9, t; R) is the same as in the previous 
Section. However, for each next generation the num- 
ber of numerical integrations increases, unless values of 
ji-i(9,t;R) are stored as a 3-dimension matrix. Thus, it 
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is convenient to find analytical solutions of the integrals 
F(ff, a) and/or G(R', a), if possible. 

2.4. Results of calculations 

Fig. 3a shows the number of photons arriving at the 
detector within an angle £ < 1°,3°, 10° per unit area 
per unit r = ct/R as a function of e = t - 1. The upper 
curves refer to the first generation, the lower - to the sec- 
ond one. We also compare here the time distributions 
obtained for the Rayleigh with those for the isotropic 
scattering. The distance detector- source equals to one 
scattering length (k = R/A- 1). 

First of all we notice that for short times (e < 0.01) the 
first generation dominates over the second one, and (as 
we can guess) over the higher ones. It can be seen from 
the formulae for the isotropic scattering (Appendix A) 

that the ratio ^p- / for any given time should be pro- 
portional to k = R/A, so that the importance of the sec- 
ond (and the higher) generation will be bigger for larger 
k. 

We can also see that the number of photons arriving 
within an opening angle £ reaches the dependence ~ £ 2 
only at later times. This reflects the fact that the ini- 
tial angular distribution of light is steep and becomes 
almost flat at times t > 1.1 or so. When comparing the 
Rayleigh curves with the isotropic ones one can see that 
the latter are slightly flatter for shorter times, as might 
be expected but become parallel to the former for later 
times. 

Next, we consider a situation when there are two scatter- 
ing processes, Rayleigh and Mie with quite different an- 
gular distributions f(a) (as discussed before). We adopt 

R - ■ ' 1 (35) 



Atot Ar Am 



and A R - A M for simplicity. 
The result for the first and the second generation de- 
pending on time is shown in Fig. 3b. There is now 
more light at earlier times than in the previous case (Fig. 
3a) due to the strong Mie scattering in the forward di- 
rections. However, the flux of the first generation de- 
creases about 3 times quicker over the considered time 
region. Although the ratio of the second to the first gen- 
eration is practically the same at e = 10~ 3 in both cases 
(< 1%). the importance of the second one is reached at 
later times when the Mie scattering is present. 

3. Point light source in the atmosphere 

Now we shall study the situation when a point light 
source flashes in a non-uniform medium, such as the at- 
mosphere. We assume that the atmosphere is composed 
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Figure 4: Vertical cross-section through the atmosphere. Lines cor- 
respond to constant values of k = R/Apo, shown by numbers, look- 
ing from the detector (at x — 0, h = 0) to a point on the line, for 
A = 18 km and A M = 15 km at the ground and the scale heights 
H R = 9 km and H M = 1.2 km. 



of two sorts of matter, molecules and aerosols, each 
having its density decreasing with height exponentially 
with a different scale heights, H R - for molecules and 
H M for aerosols, and having the corresponding mean 
path lengths for scattering at the ground Ap and Ap. It 
is not difficult to derive that the effective mean free path 
for a scattering for light travelling between two arbitrary 
points P and S equals 



(— + —) 



(36) 



hp - hs 



where h P and h s are the heights of points P and S above 
the level, where the Rayleigh and Mie scattering path 
lengths are correspondingly A 1 ^ and A^. If the source 
(point P) is at a distance R from the detector (point D 
on the ground) then the ratio k, of R to the mean free 
path length along PD equals 



k = 



R 

^PD 



(37) 



cos6 z 



H, 



VA' 



e 



-Iip/Hm 



It can be seen that increasing the distance R to infinity 
(keeping 6 Z constant) the ratio k reaches its maximum 
finite value 



kmax(6z) 



cos0 z \A R n A 



(38) 
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This situation is illustrated in Fig. 4. Here a vertical 
cross-section of the atmosphere is shown. Detector is 
at x = 0, h — and lines represent constant values of 
k corresponding to the straight path from the detector 
to the point on the line. We have adopted the follow- 
ing values: A^ = \%km,HR = 9km, A^ = 15km and 
Hm = 1.2 km. These values describe approximately the 
atmospheric conditions at the Pierre Auger Observatory 

a. 

From Fig. 4 one can also deduce that relevant values of 
k, if light sources are at distances ~ (10 — 30) km (as ex- 
tensive air showers seen by Auger) are 1/2 < k < 3/2. 
Thus, this will be the region of our interest. 




Figure 5: Geometry of the first generation in the real atmosphere. 
Light source is at P, detector at D. Scattering takes place at S . 




ground level 



Figure 6: The ellipsoids (their cross-sections are shown) show scatter- 
ing sites of first generation photons arriving at D after time ^(1 + e). 
The corresponding numbers are equal to e. Light source is at P, de- 
tector at D. 



3.1. The first generation 

As the medium is non-uniform, we cannot use the 
idea of a sphere to be crossed by the scattered pho- 
tons, as in Section 2. Now their flux will depend on 



the zenith angle 8z of the source and on the azimuth an- 
gle <p around the direction towards it. Fig. 5 shows a 
trajectory PSD of a first generation photon scattered at 
S . We want to calculate the angular distribution of the 
first generation as a function of time, for a fixed R and 

8z, trmAQ, <t>, t; R, 6z), crossing a unit area perpendicular 
to the direction towards the source. 
We notice that for a fixed arrival direction of photons 
(0, <p) and time t, the scattering point S is uniquely de- 
termined. To arrive at the detector at angles (6, <p) within 
c/Q = sin9d9d(p after time (t,t + dt), photons have to 
cross the surface da (shaded in the figure) and be scat- 
tered along a path length dx by the angle a determined 
by Eq.2. The number of such photons equals 



dn\(6,<p,i) = 



47TX 2 



dx'cosy 

— e 1 ps ■ 



(39) 



~f{a)/sn D e-^i 
As 



where x = PS , x' = S D, y is the angle between 
the normal to the surface da and the direction of the 
incident photons PS, AQo is the solid angle deter- 
mined by the unit area at D and the scattering point S 
(AQo = cosO/x' 2 ). There is no need to calculate y be- 
cause dx = x' d6 1 cosy, so that it cancels out. It can be 
shown that 



dx!_ 

~r2 



dT 



or t 2 -It cose + 1 R 
Inserting this into (40) we obtain 



ji(6,<f>,t) ■ cosO 



d 2 n 



dQdt 
f(a)cose 



(40) 



(41) 



2nA s R 2 t 2 -It cose + 1 

One can see that this formula is practically the same as 
(7) for the uniform medium, the only difference being 
in the scattering path lengths depending not only on dis- 
tances but also on the geometry. 

The height of the scattering point S necessary to calcu- 
late As , Aps and As d equals: 



= (ct- 



sinO \ 

R—\ 

sina I 



(42) 



■(cosOz cosO + sinOz sine cosip) 



We calculate numerically the time distributions of light 
dN[ ea, (t; ()/dt, arriving at the detector at angles smaller 
than ( for different zenith angles of the source. The re- 
sults, in the form of the ratio: 



Fi(t;0 = 



dN[ eal ldt 



dN l l m 



/dt 



(43) 
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are presented in Figs. 7, 8 and 9, where dN""'/dt are 
the distributions obtained in the previous section for a 
uniform medium. 















30° 




"\Royleiqh k=1/2 




' 0° 






0-01 0.1 



Figure 7: Ratio of the first generation dN" jdt in real atmosphere to 
that in uniform medium dN" m /dt as a function of time (e = ct/R - 1) 
for the Rayleigh scattering only for various values of zenith angle of 
the source. Solid lines - £ = 1°, dashed lines - f = 10°, k = Rli R PD = 
1/2. 




10 3 0.01 0.1 1 



Figure 8: As in Fig. 7 but with Mie included; k = R/A PD = 1/2. 

To compare the light flux obtained for the real atmo- 
sphere with that for a uniform medium we adopt the 
same value of k and the same distance from the source 
to detector for both cases. To understand the effect of 
a purely exponential atmosphere we start with consid- 
ering only the Rayleigh scattering (Fig. 7), neglecting 
Mie (A M = °°). Understanding the behaviour of the 
curves in this figure is easier with the help of Fig. 6. 
Each ellipse is a cross-section of a rotational ellipsoid 
with the symmetry axis determined by point D - the de- 
tector and point P - the light source. These are the focal 
points of all the ellipses. Photons arriving at the detec- 
tor at D after some fixed time t = 1 + e (in units of the 




10~ 3 0.01 0.' 1 



Figure 9: As in Fig. 8 but for two values of k. Each group of lines is 
for <T = 1°, 2°, 3°, 5°, and 10° from bottom to top (at e = 10" 3 ). 

distance PD) must have been scattered on the surface of 
an ellipsoid with the eccentricity e equal 

e = -!- (44) 

1 + £ 

The ellipses refer to e = 1(T 3 , 1(T 2 , 10 _1 and 3 • 10 _1 
keeping the right proportions. The detector field of view 
cuts only a part of the ellipsoid surface where the pho- 
tons registered after time 1 + e must have been scattered. 
We notice that all ratios F\ in Fig. 7 are smaller than 1 
and decrease with time (although those for £ = 1 ° are 
practically constant) and ratios for £ = 10° are larger 
than those for f = 1°. All this becomes clear when in- 
specting Fig. 6 and the corresponding scattering sites. 
For example - the constancy of F\(t; £ = 1°) for times 
t - 1 = e = 1(T 3 -r 1CT 1 reflects the fact that the scat- 
tering takes place very close to point P during all this 
time and starts to move away from it (higher in the at- 
mosphere) only for larger times i.e r > 1.3. Since the 
scattering path length in the uniform medium is chosen 
equal to the effective path length in the atmosphere Apd 
(Eq. 36), we have that A P > A PD and the scattering prob- 
ability at P is smaller in the exponential atmosphere. 
Let us consider now a more realistic atmosphere with 
both processes, Rayleigh and Mie at work. Fig. 8 shows 
time dependence of the ratio F\(t; £) for £ = 1° and 10° 
and zenith angles 9 Z = 10° 4- 75°. Let us take a closer 
look at the case £ = 1° and 6*z = 75° (upper solid curve). 
It may seem strange that the ratio F\ increases since 
in the case of the Rayleigh scattering only it decreases, 
although very slowly. We have checked that a similar 
slow decrease takes place if the scattering is only Mie. 
The behaviour of the curves in this figure would be hard 
to understand without the presentation of the scattering 
sites in Fig. 6. It can be seen that for e smaller than a 
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few xlO 2 the scattering sites are close to the source at 
P. Thus, we may approximate the ratio F\ as follows 



Fi 



fR,fM d 



I*. 

'Ion 



J M 
) M 



i , 2_p_ Jm_ 

lfi 7 

1 , ^PD J M 
1 "•" iM 7 

A PD /ft 



(45) 



where / K(M) are some effective angular distributions of 
photons scattered by Rayleigh (Mie) on the cut surface . 
As time increases, none of the A's changes much. How- 
ever, the typical scattering angles increase, what affects 
much more f M than f R , so that f M / f R decreases. Since 



A P A PD 
A P A PD 



(46) 



the numerator decreases by a smaller factor than the de- 
nominator so that the ratio F\ increases. It can be seen 
from Fig. 6 that for e > 2 • 10~ 2 the scattering angles 
of the registered photons do not change much (the de- 
nominator stays constant) but now Ap and A R have to 
be substituted by /I? and Af, where S is an effective 
scattering point with growing height. As H M < H R , the 
ratio /I* /Af decreases and so does F\ . 
A different behaviour of F\ (t; f = 10°) can be explained 
also with the help of Fig. 6. At first F\ decreases (the 
smaller 6 Z - the stronger decrease) because the detec- 
tor field of view cuts out a growing part of the deep 
atmosphere where the scattering is strong in the real at- 
mosphere. At e > 0.02 F\ starts to increase for the 
same reason as just described in the case £ = 1°. It 
must finally decrease since the scattering takes places 
further and further behind the source, where A R and A M 
are growing in the real atmosphere. 
In Fig. 9 we show F\{t; for 9z = 60° for several inter- 
mediate values of and for two values k = 1/2 and 1. 
Note that changing k = R/Apo must result in changing ft 
- the distance to the source. This is the main reason why 
the curves for k = 1 are lower than those for k = 1/2. 
It can be seen from Fig. 5 that for 6z = 60° (the scales 
on both axes are the same so in the figure the angles are 
correctly represented) the distance ft is much shorter for 
k = 1/2 than for k = 1, implying that the correspond- 
ing heights of the source differ considerably. Inspecting 
Fig. 6 we can estimate that for e 0.01 and £ = 1° the 
curves for k = 1 should be down with respect to those 
for k — 1 /2 by a factor 



^■P,D i A PgiD 



(47) 



where P1/2 and P\ are positions of the source referring 
to k = 1/2 and 1 respectively, with the meaning of all 



A's as defined before. In our example A PlD ^ 33 km, 
Ap x 113 km, Ap i/2 o — 15.3 km, Ap i/2 25.6 km so that 
the above factor equals 0.49 whereas the exact ratio 
from Fig. 9 - 0.48. 

3.2. The second generation 




Figure 10: Ratio of the second dNildt to the first generation dN\ Idt 
as a function of time (e = ct/R - 1) from a flash (at t = 0) of a point 
source at zenith angle 9z = 75° in the real atmosphere. Fluxes are 
integrated within f = 1° (solid lines) and £ = 10° (dashed lines). 
Three curves for each f refer to k = 1/2, 1,3/2 (from bottom to top). 
Dotted line refers to a uniform medium with Rayleigh scattering only, 
forfc= 1,£= 1°. 

We proceed like in the case of the first generation 
(Fig. 5), but now point S refers to the second scatter- 
ing. Photons scattered only once arrive at the surface da 
from all directions according to ji(0i,4>i,h;x), where 
t\ = t — x' /c. Thus, the number of photons dn2 inci- 
dent on da and scattered for the second time towards 
the detector (to arrive there within dQ.o{Q, <p) after time 
(f, t + di) equals 

dri2(6, <p, t; x') = I ji{6i,(f>i,t;x)dQ.\dtdacosy- 



dl 

• —f(a 2 )dQ. D e *sd 
As 



(48) 



where the integration has to be done over full solid angle 
(0 < 0i < 2tt, < 61 < n), dl is the path length for the 
second scattering to occur (cosydl = x 1 dff) and a 2 is 
the angle of the second scattering. Now the pair of vari- 
ables, 6 and t, does not determine uniquely the positions 
of the second scattering, since the times t \ elapsed from 
photon emission to their arrival at S (or strictly speak- 
ing, at da) have some distribution. However, t] can not 
be smaller than x/c, thus x' max = ct — x. Expressing x' nmx 
as a function of 6, t and ft only we obtain 

{ctf - ft 2 



X — 



2(ct-Rcos9) 



(49) 
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and the distribution of the second generation equals 



d 2 n 2 (6, (p, t) 



dQ dt 

r 

I d9\ j\ sinOi 
Jo 



I dx I ( 

Jo Jo 



dtp- 



e *sd cost! 



(50) 



Finally, the number of photons within an angle £ 
equals 



dN 2 (t;Q 
dt 



r*2jr r*7T 

Jo Jo 



d 2 ti2 
dQdt 



sinG d6 dip 



(51) 



The above integrals have been calculated numerically. 
The ratio / ^ as a function of t - 1 in the real atmo- 
sphere, where the source is at Gz = 60° at the distance 
R - 34 km corresponding to k — 1 , within two open- 
ing angles of the detector £ = 1° and 10°, is presented 
in Fig. 10. For comparison we have also drawn there 
the same ratio for the uniform medium, for k — 1 with 
Rayleigh scattering active only (following from Fig. 3a 
for f = 1°). First we see that the contribution of the 
second generation increases with time as it should be 
expected. Next, it is considerably larger for the uniform 
medium than for the real atmosphere, which may not 
be so obvious at first sight. In fact, at e = 10~ 2 -f 10" 1 
the curve for the uniform medium is ~ 2.3 times higher 
than that for the real atmosphere. This factor equals ( 
for simplicity we neglect the operators Jj) 



n = 



N u 2 ni {R) N r 2 eal {R + M) 
N l l ni {R)' ' N[ eal (R + M) ~ 
N[ eal {R + M) N r 2 ea '(R + M) 



(52) 



N" ni (R) 



N" ni (R) 



The numerator ~ = 0.29 as explained in the 
previous paragraph. The denominator refers to photons 
scattered exactly two times. For small e the scatterings 
must take place close to the source or along the field of 
view = 1° in our example). Thus we should have that 



N r 2 eal (R + M) i A 



N l 2 lm {R) 



084 



(53) 



Our calculations show that 0.13 for e = few 

Af 

xlO^ 2 so that j] ~ =i 2.2, in agreement with the 
exact calculations. 



4. A moving point source - a shower 

So far we have been interested in the light signals 
from a flash of a stationary point source in a detector 



at a given distance, within a given field of view, as a 
function of time. 

Now, we shall consider a moving point emitting light 
on its way through the atmosphere. This is a reason- 
able model of a distant cosmic ray shower exciting at- 
mosphere what results in emitting isotropic fluorescence 
light. Moreover, the Cherenkov light produced by ~ 1/3 
of shower electrons J5t], propagating alongside the par- 
ticles will be scattered to the sides by the Rayleigh and 
Mie processes in an anisotropic way. Thus, in principle, 
a shower is a moving source, emitting light anisotrop- 
ically. Our aim is to find an instantaneous optical im- 
age of the shower, produced by the scattered light at the 
same time as the shower image in the direct (not scat- 
tered) light. 

The image obtained for a small integration time Af is a 
radially symmetric light spot, corresponding (roughly) 
to the emission distribution at the shower lateral cross- 



section Ill2ll . It is produced mainly by the direct light. 
But at the same time as the direct light, there arrive also 
photons produced by shower particles at an earlier stage 
of shower evolution and scattered in the atmosphere. 
As shower parameters follow straightforwardly from the 
amount of the direct light, it is important to calculate 
this effect and to take it into account when deriving the 
flux of the direct light. 

Let us assume that we know the geometry and the tim- 
ing of the shower. If the source is at point Q at time t 
and emits C photons per unit length, the telescope will 
record the following number of the direct photons at 
time t + R/c, in the time bin Af: 



cAt A s_ 

An = C - -Jp^)^ QD 

1 - co so R z 



(54) 



where 6 is the angle between the shower direction and 
the direction from point Q to the detector D, Aqq is the 
effective mean free path for attenuation along the dis- 
tance R, A is the diaphragm area determining the collec- 
tion solid angle and fp(6) is the angular distribution of 
the emitted light. 

In the case of isotropic fluorescence light fp{5) = l/4n 
and C equals to the number of photons produced per 
unit length (~ 4mT x , per one electron). In the case of 
the scattered Cherenkov photons, or those from a laser, 
fp(6) = j4p(l + cos 2 S) for the Rayleigh scattering (or 
a corresponding distribution for Mie) and C = N/Aq, 
where Aq is the mean free path (in length units) for 
the process in consideration at point Q and is the 
number of Cherenkov photons at the observed shower 
level propagating (approximately) in the direction of the 
shower. The number of the scattered Cherenkov pho- 
tons to the fluorescence ones is typically 10 - 15%, and 
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a lot of the former are scattered by Rayleigh, what is 
not far from isotropic. Thus, with a good approxima- 
tion one can treat a shower as a moving isotropic light 
source. This is what has been adopted in this paper: 

Simultaneously with the direct light some photons pro- 
duced earlier (above point Q) will also arrive at the tele- 
scope. However, as we already mentioned, they must 
have been scattered on their way to the telescope. Our 
earlier calculations [10] showed that if the distance to 
the shower is not very much longer than the mean scat- 
tering path, then the main contribution to the scattered 
light is due to the photons scattered only once (the first 
generation). Here, we will calculate the fraction of the 
total light consisting of the first and the second genera- 
tions. 

4.1. The first generation 

Let us consider photons produced at point P, along a 
shower path element dl (Fig. 11). The first generation 
photons arriving at the detector D at the same time as 
the direct photons produced at Q must have been scat- 
tered on the surface of the rotational ellipsoid with the 
focal points at D and P, with the DP line being the ro- 
tational symmetry axis of the ellipsoid. 
As we have already mentioned, the main interest of the 
shower experiments is to measure the direct light. Thus, 
we are interested in the light arriving from the direction 
of point Q and the region around it. In practice (as it is 
in The Pierre Augere Experiment) the light signal arriv- 
ing simultaneously with the direct light from point Q is 
measured within a certain angle f around the line DQ. 
The cone with the opening angle £ cuts out on the sur- 
face of the ellipsoid a region of our interest (the shaded 
surface). It is at this part of the ellipsoid that the pho- 
tons produced at P have to be scattered to arrive simul- 
taneously with the direct light produced at Q, within an 
angle smaller than f with respect to the latter. Thus, to 
find the contribution of the first generation to the direct 
light one has to integrate the number of the scattered 
photons over the cutout surface, and then integrate the 
result over the distance QP = I (the upper part of the 
shower). 

The proportions in Fig. 1 1 have not been preserved. 
For example the actual distance from point P to the 
detector at D is much longer than / = QP (although 
in principle the integration over / should go far up the 
shower, the distances contributing to the total scattered 
light are rather close to point Q). Also the angle £ does 
not need to be larger than a few degrees, but on the fig- 
ure it is much larger for the sake of clarity. 



shower 




Figure 1 1 : First generation from a shower. The shaded surface shows 
sites of scattering of photons produced at P, scattered once and arriv- 
ing at detector D within its field of view simultaneously with direct 
photons produced at Q. PQ+QD = PP\ + P t D. 



Let us first find the contribution to the first generation 
from photons produced by the shower at point P along a 
path length element dl, assuming that P lies outside the 
detector field of view. 



a)0 P >f(Fig. 11). 

(dni\ cCdl r 9 ^ 



f(a) sinO cosO 



t' 2 - It' cosQ + 1 



Jo Ap, 

l+R 



-d(f> 



(55) 



where t' = tg(a/2) 



t' —cos9 




Figure 12: Fluorescence light produced at point P and scattered on 
the shaded surface (inside field of view of detector D) of the ellipsoid 
arrives at D simultaneously with the direct light produced at Q. 
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the direction to the source at P), and R' = PD = 
V/? 2 + I 2 - IRlcosd. 

The azimuth angle (p is measured in the plane perpen- 
dicular to the axis PD and = refer to points on 
the shower-detector plane (Fig. 1 1 plane). The length 
PP X can be found from the triangle PP X D: PP X = 
R'sinO/sina and PiD = l + R-PP v 

b) 9 P < (. (Fig. 12). 

It is clear from the figure that now the integration limits 
of 6 and (p are different so that one has to add to the 
integral in case a) the following term 



cCdl 
2nR' 2 



J f(a)sin6cos6 
d0\ 



1 i 

Jo 



t' 2 - It' cosd + 1 



(56) 



-i 1 + 1 ) 



and the lower limit of in the term a) changes sign, so 
that Op - £ -> C ~ Op. To find the total flux of the first 
generation arriving at the same time as the direct light 
from point Q one has to integrate over / the contribu- 
tion from case b) from to l max = Rsin£/sin(S - f), 
(where 6 > f), and add to it the contribution from a) 
integrated from l max to oo. (It turns out that contribution 
from points with / 0.1R are negligible). 
If there are two scattering processes the first generation 
means the sum of the number of photons scattered by 
Rayleigh and those by Mie and then ^ 



A,; 



f(a) f M (a) 



and in the attenuation exponent A PPl and A PiD are the ef- 
fective mean free paths for both processes. 

4.2. Analytical calculation of —jp- 

By making some approximations we have found an 
analytical solution of the shower image in the first gen- 
eration light scattered by the Rayleigh process. 
We shall treat separately light produced within the de- 
tector field of view ^ and that outside it, —^ L , starting 
with the former (Fig. 12). As the scattering points lie 
close to point Q we assume that Ap { = Aq. For the same 
reason we assume that the exponential factor describ- 

R 

ing light attenuation equals e A o° . Our main assump- 
tion, however, consists in integrating over the symmet- 
ric (with respect to axis PD) part of the ellipsoid so that 
we could solve the integrals over (/>, 6 and finally I. Ex- 
pressing scattering angle a as a function of 6 and r' (Eq. 



2) we obtain that the contribution from a shower path 
element dl equals 



V dt / R' 2 8tt 



(57) 



r """ , i 2sin 2 2sin 4 6, sinO cosO d9 

Jo ' y y' 2 ' y 

where y' = t' 2 - It'cosQ + 1, B = ^e' 7 ^ . 

We adopt 9 max = £ being the mean of the limiting values 

of 6 in the exact integration. 

In the air shower experiments it is the direct light which 
should be measured so that the viewing angle f is small 
(a few degrees) and so are angles 6. This leads to t'- 1 = 
e <K 1. Taking all this into account the result of the 
integration over is 



\~~dT) ~ 



ln- 



cCe a qd 3 
Sn' 

1 



(58) 



+ 1 



(x 2 + l) 2 



dl 
^ 2 



where x = e/g. From geometrical considerations it can 
be derived that 



dl_ 

R 12 



1 



ecosS 



1 - cosS (1 - cosS) 2 
dr 1 £dx 



dY - 



(59) 



R{\-cos8) R{\-cos6) 
where we have dropped the second term ~ e. 

dn R 

To obtain it remains to integrate (59) over x (includ- 
ing (60) from to x max , corresponding to point P lying 
on the edge of the viewing cone (9p = f). Again, as- 
suming that £ «: 1, it can be derived that x max tg^f 
(for 6 > £ ). The result of the integration is the following 



dt AuRAqH - cos6) 
3 



(60) 



S' 8' 
sin§ - 8tg — • In sin 
2 



where 6' = S - f . 

To calculate the contribution from the outer part of the 
shower i.e. outside the field of view we keep the previ- 
ous approximate assumptions. However, since now the 
scattering surface looks differently (Fig. 11) we make 
another simplifications. We assume that all scattering 
angles equal 6 so that f(a) = f(S) = const. Next, as £ 
is small we assume that all 8 = 6 P . Integration over Q fl , 
the solid angle of the field of view, is then reduced to 
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Figure 13: Second generation from a shower. Photon, emitted at P 
must be scattered for the second time (point Po) within detector D 
field of view. First scattering (Pi ) must occur on the surface of the 
ellipsoid with focal points at P and Pi. Direct light from Q arrives at 
D at the same time: PQ+QD = PPi +P1P2 + PtX>. 



multiplying the integrand for 6 — 6p by Q.o = ji( 2 . Thus 
we have 

J dn Li \_cCe~^dl C f(a)dQ(0, d>) cosd ^ 
\ dt I ~ 2jtAqR' 2 ' J nD t' 2 - 2t'cos6 + 1 ~ 



cCe l w> f{5) ntfcosOp dl 



(61) 



2tt A q t' 2 - It'cosB + 1 R' 2 

It can be derived that the denominator in the above ex- 
pression equals 

R' 2 (t' 2 - It'cosB +1)^(1+ sin 2 5 - cos5)l 2 (62) 

Assuming^] that cosOp — 1, we obtain 



(d< ut \ 
\ dt I 



cCe 'eo f R (5) 



In 



Since /, 



1 + sin 2 5 — cosd 
~ R i 



r°° di 



-» we have 



mm sin(S-Q 
R ^•f~'^~ inn r 1 (\ _i_ ™„2.$ 



dn K ul cCe A Q° £ 3 (1 + cos 2 5)sin{6 - 
dt 2nRAn 16 1 + sin 2 5 - cosd 



(63) 



(64) 



1 Taking cosOp = 1 - 4 sin5( ^ ) 2 needs to allow for light attenuation 
along / in the integration over /. The correction to dn ou ,ldt is a few 
percent, being negligible for the total dn'J/dt. 



dn R 



dnj | dn R Ui 
dt dt dt 



(65) 



Expressing the number of photons An* = ^-Af as the 
ratio to that of the direct light (Eq. 54 for f(S) = 1 I An) 
we finally obtain 



A "i „ 3 R 

A«() 16 Aq 



5' 6' 
36 - sind - &tg— • In(sin-) 



+ 2 



sin5'(\ + cos 2 6)(\ — cosd) 



1 + sin 2 5 — cosd 



(66) 



where 5' — 6 - £ (in radians). Terms oc £ 2 and of higher 
order have been neglected. We can see that, with the ap- 
proximations adopted, the ratio An*/A«o is proportional 
to R(/Aq reflecting the fact (as we shall see below) that 
most scatterings take place close to point Q since R( 
is proportional to the shower path segment seen by the 
detector. For £ «: 5, 5' 6 and the dependence on 
6 separates from that on other parameters. It becomes 
also obvious that it must be different for different an- 
gular distributions /(or), of photons at scattering. As 
f R (a) oc 1 + cos 2 a is quite different from f M {a) which 
is peaked in forward directions, we expect another de- 
pendence on 5 for the latter. 

While deriving the analytical formula (66) we have also 
assumed that the length of the shower segment cut out 
by the detector field of view is small when compared to 
the distance R. For £ <k 1 it is fulfilled for almost all 5, 
apart from the case when 5 is close to or n. We have 
calculated analytically the ratio A«*/A«o for 6 = n - y 
for y — > 0. However, as the derivation is lengthy and 
(as it will be seen later) agrees with the exact, numeri- 
cal results only at 6 = n, we present here only the final 
result: 

An? 9 / 7 N 

An J^ n 16 v 3n ' 
These values are marked as stars (for 5 = n) in Fig. 15. 
We have not tried to calculate Anf /dt because the func- 
tion f M {a) is not well known, depending on the sizes 
of the aerosol particles and, anyway, our analytical ap- 
proach is approximate. Thus, it seems better to calculate 
it exactly numerically and, being led by our analytical 
solutions for Rayleigh, find an appropriate parametrisa- 
tion of the numerical results. 

4.3. Numerical (exact) calculations of A«i/A«o and 
A« 2 /A«o 

The ratio of the first generation to the direct light ar- 
riving simultaneously is obtained by numerical integra- 
tion of point contributions over I (shower track above 
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Figure 14: Ratio of the first generation to the direct light as a function of distance R to point Q on the shower in units of the scattering length Aq 
at that point, a). Rayleigh only b). Mie only. Curves correspond to ( = 1°, 3°, 5° (from bottom to top) and 6 = 90°. Various point shapes refer to 
different distances R=\2, 24, 32 km. Lines are power law fits. 



point Q ) and dividing the result by A«o- Fig- 14 shows 
the ratio for the Rayleigh a) and Mie b) scattering as a 
function of k Q = R/A Q for f = 1°,3° and 5° and for 
several different distances R. It can be seen that for the 
Rayleigh case the ratio is proportional to k Q and f , as it 
has been derived analytically. For Mie the dependence 
follows a power law, with the indices depending slightly 
on However, in each case there is practically no de- 
pendence on the distance R itself. These results are for 
5 = 90°. 

The dependence on 6 is shown in Fig. 15a,b for fcg = 
1 and £ = 1°, 3° and 5°. One can see the difference of 
the behaviour of this dependence between Rayleigh and 
Mie. In the figure a comparison of our analytical calcu- 
lations with the exact ones is also shown. The biggest 
difference reaches some 10% at large 6. There is no 
normalisation there. We find this agreement quite satis- 
factory. 

Nevertheless, as the agreement is not perfect we have 
parametrised the exact numerical results. For the 
Rayleigh scattering the ratio A«*/A«o can be expressed 



as 



Anf _ 
— 1 = 0.024 • k Q ■ (g R (x) 
An 



(68) 



where x = 5/100°, f is in degrees, and 

g R (x) = 0.112 + 1.86* - 1.33x 2 + 0.383x 3 (69) 

For Mie the factorisation is not as complete as for 
Rayleigh, and our fit to the numerical results is the fol- 
lowing: 



Anf 
Atto 



(70) 



kQ x0.93-0.04f / f x 0.93+0.05, 



= 0.096 •(-%) •[£-) 



g M (x) 



x < 0.6 g M (x) = x(7.76x z - 11.25* + 5.625) 



x > 0.6 g M (x) = 1 



(71) 



These fits are also shown in Fig. 14 and 15. It is seen 
that they do not deviate from the exact values by more 
than a few percent. 

The second generation has been calculated numerically. 
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Figure 15: Ratio of the first generation to the direct light as function of angle S between shower and line of sight, a). Rayleigh only b). Mie only. 
Solid lines - exact numerical calculations, dotted lines and stars in a) - analytical (approximate) calculations, dashed lines - our parametrisation of 
the numerical curves. All curves are for Icq = R/Aq = 1. 



Fig. 13 should be helpful for understanding the calcula- 
tion. The thin solid line shows a path of a photon pro- 
duced at P, scattered at points Pi and P2 and arriving at 
D. The point P2 of the second scattering must lie some- 
where inside the viewing cone. Fixing it at a distance rj, 
from the detector corresponds to fixing the value of the 
remaining photon path length PP1P2. We have that 

PPi + PiP 2 =R + l-r 3 (72) 

where I = PQ and R = DQ. 

In order to assure that the photons scattered twice ar- 
rive at D at the same time as the direct light emitted 
at Q, point Pi must lie anywhere on the surface of the 
ellipsoid with focal points at P and P2. Thus, to calcu- 
late the number of photons produced at P and scattered 
twice one has to integrate the contributions like in (49) 
over the position of the first scattering (the whole sur- 
face of the ellipsoid) and then integrate the result over 
the volume of the viewing cone. Finally, the integration 
has to be performed over the distance I = PQ. 
The ratio A«2/A«o, including Rayleigh and Mie scatter- 
ings as a function of k Q = R/Aq for 6 = 90° and as a 



function of 6 for k Q = 1 is presented in Fig. 16. Various 
point signs refer to different distances R. As in the case 
of the first generation the ratios depend on R/Aq rather 
than R itself. However, now they are proportional to the 
second power of it. We have parametrise the obtained 
numerical results as follows 

= 1.2- 10->W"Veg) (73) 
A« u 90 

5. Discussion of shower results 

5.1. Dependence of A "^" 2 on height 

As we see the numbers of the scattered photons, A«i 
and A«2, depend on the angle £ within which they ar- 
rive at the detector. This angle has to be chosen in such 
a way as to encompass the total direct signal from the 
shower. It is obvious that for closer showers f has to be 
larger and for those more distant - smaller, so that R£ 
- meaning the lateral spread of the direct light - should 
remain constant. Our approximations show that in this 
case (R£ = const) nf/no depends only on height (by Aq) 
and 6. The power indices in the dependence of Am /A«o 
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Figure 16: a). Ratio of the second generation to the direct light as function of Icq. b). Ratio of the second generation to the direct light as function 
of angle 6 between shower and line of sight. Points - results of numerical calculations, lines - our parametrisation. Different point signs refer to 
different distances (R = 12, 24, 32 km); kg = R/Aq = 1. 



on R and £ are close to 1, whereas A«2/A«o ~ (R • <T) 2 
(roughly), so that their dependence on the distance R it- 
self should be weak. Fig. 17 presents the contributions 
of the individual components of the scattered light as a 
function of height, for R( = 30 km- deg, 6 = 90° and two 
(quite different) values of R: 12 and 36 km. The depen- 
dence of the sum "^" 2 on R is very weak, indeed, the 
more so as the ratios Anf /Ana and A«2/A«o change in 
opposite directions. The chosen value R( = 30 km deg 
corresponds to a shower lateral radius of ~ 500 m, equal 
to ~ 5 tm (Moliere radii). This may seem too large to 
contain the total fluorescence signal. However, deep in 
the atmosphere it is the shower Cherenkov light what 
determines its lateral dimensions and it goes up to 5 ym 

a. 

As the typical viewing angles d are not far from 90° one 
can draw a conclusion that the maximum contribution 
of the MS light is ~ 14% for RC, = 30 km deg and the 
adopted atmospheric parameters. 



5.2. Influence on a reconstruction of shower parame- 
ters, Eq and X max 
Let us estimate how not allowing for the MS effect 
would affect a reconstructed value of a shower primary 
energy Eq. The total signal would then be treated as 
the direct light, therefore the primary energy would be 
overestimated. To estimate how much it would be let us 
consider a typical (for Auger) shower with Eq = 10 19 
eV, falling to the ground 25 km away from the light de- 
tector, at 6z = 40° and azimuthal angle 90° (angle be- 
tween the direction from the detector to the shower core 
and the projection of the shower axis on the ground). If 
the primary is a proton the mean depth of its maximum 
equals X max 790 g cmr 2 , corresponding to the height 
above Auger level (~ 860 gem -2 ) h max =t 3.2km. From 
Fig. 17 we see that the contribution of the scattered 
light is ~ 4.5% if R( = 30 km deg, what corresponds 
to £ = 1.2°. However, in the Auger telescopes the indi- 
vidual pixel has a diameter of 1.5°, so that it would be 
necessary to choose also the neighbouring pixels around 
the one with the stronger signal to collect the total direct 
light. Thus, we choose f = § • 1.5° = 2.25°. Within this 
angle the MS light grows to ~ 9%, overestimating the 
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Figure 17: Ratio of the first two generations to the direct light as a 
function of height of the observed shower. Solid lines - R = 12 km, 
dashed lines R = 36 km. Curves denoted by "Rayleigh" and "Mie" 
refer to the first generation, by "n2" - to the second one. If a constant 
lateral dimension is chosen (R( = const, S = 90°) the ratios do not 
(practically) depend on shower distance R. 

number of particles at X max by this value. For shower 
parts near the ground, where shower age equals 1.25, 
this number grows to ~ 27%. To find the shower pri- 
mary energy Eq (by integrating the energy deposit over 
depth Xq) it is necessary to extrapolate the shower for 
larger ages, so that the overestimation from these parts 
would be even larger. We estimate that not allowing for 
the MS effect would increase the reconstructed £o by 
~ 15%. 

Not allowing for the MS light would also lead to an 
overestimation of the depth of the shower maximum 
Xmax. Fitting a new cascade curve to that from the pre- 
vious example with the increased values according to 
the MS contribution, we obtain that the overestimation 
equals tSX„ mx — 35 g cm~ 2 . With the difference between 
the depths of an average proton and iron showers be- 
ing ~ 100 gcm~ 2 , this is not a negligible number for 
drawing conclusions about the mass composition of the 
highest energy cosmic rays. 

5. 3. Contribution of higher generations 

At first sight it may seem odd that it is only the first 
two generations that contribute to the scattered light. 
One might expect that the main contribution should 
come from the generation with the number closest to 
the value R/Aqd, according to the Poisson distribu- 
tion (although we have shown that the first two gener- 
ations scale with R/Aq rather than with R/Aqd). For 
Ar = \%km and Am = 15 km and a shower close to 
the ground we have Ago = %.2km . Thus, if the dis- 
tance R = 24 km then R/Aqd — 3, so should one expect 
the third and the fourth generations to contribute even 



more than the first one? The answer is no, since the 
scattered photons should fall into the detector field of 
view, so that it is not the path length of a scattering by 
any angle (what Aqd means ) that counts but by some 
particular values of it. This decreases the probability 
of such a scattering and can be viewed as an increase 
of the (actually) effective mean free path length A e ff for 
the process we are interested in (i.e. photons arriving 
at the detector). One could estimate its value from the 
ratio ni/«o or «2/"o adopting for this aim a Poisson dis- 
tribution of the number of contributing generations. The 
two obtained values differ typically by 30-60%, so that 
an estimation of the share of the remaining generations 
on the basis of this uncertain A e ff is uncertain as well. 
However, let us take an example: R — 25 km, eleva- 
tion angle 15°, £ = 2°, what are typical conditions for 
registered showers in the Auger Observatory. The two 
values of the mean number of (effective) scatterings are 
0.034 (from «i/«o) an d 0.061 (from «2/"o)- Assuming 
the larger one, one gets 2.7% for a fraction of the third 
and higher generations. 

However, as our earlier work for a uniform medium 
ifioll b) shows, where the third generation was calculated, 
its share was two times smaller than that deduced from 
the Poisson model and the ratios «i/no and «2/»o> so 
that this estimation may be too high. Thus, we estimate 
that the contribution of the higher number generations 
should be small, at least for the conditions of Auger. 

5.4. Comparison with other work 
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Figure 18: Ratio of the scattered to direct light An/Ang integrated 
within angle f , for shower 24 km away from the detector, at elevation 
angle 10°. Our results (first and second generation) are represented 
by thick solid lines for 6 = 90° and thin solid lines - for S = 30° and 
150° (second generation). Parametrisation from Monte-Carlo works: 
Roberts - dashed line and Pekala et al - dotted line. 

The problem of the influence of the scattered light on 
the shower image in the direct light was first undertaken 
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by Roberts JjJ. He treated it by Monte-Carlo simula- 
tions of tracks of individual photons emitted by verti- 
cal showers. Such shower (although, strictly speaking, 
non-existing) enabled the author to introduce a detector 
in a form of a ring around the shower as its symmetry 
axis, increasing dramatically the number of registered 
photons. Although the author appreciate that the singly 
scattered photons dominated the total MS signal, his 
simulations went to more than five scatters, including 
both Rayleigh and Mie. He obtained a parametrisation 
of the ratio of the MS signal to the total signal (MS + 
direct) in the form 



K = 77 A(OD xa-R 



1/2 



^1. 1^0.68 



(74) 



where OD is the total optical depth (molecular and 
aerosol) between the detector and the light source, a - 
total scattering coefficient (in mr l ) at the source. In our 
terms OD - R/A 0D , a = y- = 4- + i, where A 0D 

Aq A g A g 

is determined by formula (36) and Aq depends only on 
the height in the stable atmosphere, so that (74) would 
mean the following 



/ 1 N 0.68 
K-R.f'l—LJ) 
\AqdAq) 



(75) 



This parametrisation differs from those obtained in the 
present work although not very much. It has a weaker 
dependence on the angle £, which in our case is as 
(about) ~ f 1 . The smaller power index at A^ M (~ Aq 1 
in our case) is to some extent compensated by the term 

1-0.68 
QD ■ 

Another work, dealing with the present problem was 
that of Pekala et al 111 ill . The calculation method was 
Monte-Carlo as well. This paper did consider inclined 
showers so that the cylindrical symmetry, as in the pre- 
vious case did not apply. Instead large "packages" of 
photons were assumed to behave in the same way, ac- 
cording to simulated points and directions of their scat- 
tering. After extensive simulations the authors arrived 
at a parametrisation of the ratio M of the MS signal to 
the direct one as follows 



M oc 



R 

QD 



A 



(76) 



with G = 5.43 km. For a one component exponential 
atmosphere with the scale height G we would have 



1 

To 



(77) 



Indeed, it is not difficult to check that an effective scale 
height of the two-component atmosphere used by the 



authors for 2 < h < 6 km is about 5 km what approxi- 
mately equals G. Thus, their parametrisation would be 

as 



M=oc R( 



1 



AqdAq 



1 



(78) 



It is not far from our derivation for the first generation, 
where the term ^ does not appear. A comparison of 
the ratio of the total MS signal to the direct one as a 
function of £, obtained in this study with that of Roberts 
and Pekala et al is presented in Fig. 18. The agreement 
is quite good for £ < 4° but for larger ( the other authors 
seem to have underestimated a little the effect. 



6. Summary and conclusions 

Previous studies of the influence of the scattering of 
photons on their way from the shower to the detector on 
the shower image used Monte-Carlo simulations. The 
effect has been called "the multiple scattering" since 
it was not clear how many times the photons arriving 
simultaneously with the direct light actually scatter. 
Although Roberts J9[] appreciated the importance of 
photons scattered one time only he went on with 
simulating till several scatterings. 
The main idea of this study was to divide the scattered 
light into separate generations: the first generation 
being photons scattered exactly one time, the second 
one two times, and so on. This enabled us to find 
analytically the angular distributions of photons as a 
function of time for the first generation which turned 
out to be the most important one. We also showed that 
any next generation can be found by some numerical 
integrations of the previous one. These results were 
obtained for a point source in a uniform medium 
as well as in the actual atmosphere (although more 
numerical integrations were needed in the latter case). 
By considering separately the first generation it was 
easy to understand the differences between the two 
cases. The ratio of the second to the first generation 
is very small for small times after the light emission 
which will turn out to be relevant when considering 
the scattered light from extensive air showers what is a 
final aim of this study . 

In particular, this aim was to calculate the angular 
distributions of scattered photons arriving simultane- 
ously with the direct (not scattered) light emitted by 
an extensive air shower and find relevant parameters 
determining these distributions. The solutions for a 
point source were useful here since a shower could be 
treated as a series of consecutive (in space and time) 
point sources. For the first and the second generations 
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(the most important) we have explained where actually 
the scattering points took place, so that we were able to 
integrate over all of them. 

We have derived an analytical expression (although 
approximate) for the angular distribution of the first 
generation photons, scattered by the Rayleigh process, 
fitting quite well the exact numerical distributions. The 
analytical solutions have led us to a proper choice of the 
independent variables for a final parametrisation of the 
numerical results. For example the number of photons 
in the n-th generation turns out to be proportional to 
(R/Aq)", at least for n < 3 (although for Mie with some 
correction). Knowing the positions of the scatterings 
we were able to show that there should be a dependence 
on the viewing angle 6, not taken into account in the 
considerations of the previous authors. Moreover, 
this dependence should be different for the Rayleigh 
and the Mie scatterings. Thus, we parametrised the 
first generation separately for these two processes. 
Concerning the second generation we combined all 
four possibilities of the two scatterings (RR, RM, MR 
and MM) to give just one parametrisation since its 
contribution was rather small. 

From our parametrisation it follows that for R£ = const 
(corresponding to a constant observed lateral dimension 
of a shower) the contribution of the MS component to 
the direct light depends practically only on the height of 
the observed shower element and on the viewing angle 
(the latter dependence being rather mild). A typical 
overestimation of the primary shower energy (if the 
MS effect was ignored) would be ~ 15% and of the 
depth of shower maximum tsX max =t 35 gcirT 2 . Our 
analytical and numerical approach gives similar results 
as the previous Monte-Carlo work of other authors, 
although the parametrisation obtained in this paper is 
different. Only for larger opening angles > 4°) the 
prediction of the MS contribution obtained here is a 
bit larger, what would have some implications for a 
reconstruction of very close showers. 
Summarising we conclude that the MS effect is not 
to be neglected while reconstructing parameters of 
the extensive air showers. To be properly allowed for 
a good knowledge of the atmospheric conditions is 
necessary [13] for determining values of Aq, as the 
shower traverses the atmosphere. 
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Appendix A. Isotropic scattering of light in a uni- 
form medium 

It may be of some value to know how much the fluxes 
j R (0, t; R) for the Rayleigh scattering differ from those 
when the scattering function f(a) = const, so when the 
scattering is isotropic. To calculate it we proceed as in 
Section 2 putting f(a) = j- in Eq. 6. We obtain 



„-kr 



A s (e,t;R) = 



8n 2 AR 2 t 2 -2tcos0+1 



(A.l) 



and 



dt 



c e 



-kT 



(r 2 + 



\6nAR 2 T 2 
1 - 2tcos( + 1 
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)-2t(1 -cosQ 



The second generation can be obtained by proceeding 
as in Section 2.2, giving the result 



where r' 



Cl-X 

R' 



R' 



yx 



a 



IRcosOx' + R 2 and 



x' = ct-R' = f 1 -„. To find ^f^- one must 

max 2 T-cos9 dl 

integrate the flux over 6 numerically. How does the ra- 
tio of the two generations depend on the distance ? We 
can represent the first generation (A.2) in the form 



<* -7/1 



dt 



R 3~ 



(A.4) 



20 



where f\ depends on t and £ only. 

The flux of the second generation equals 



-kr 



(A.5) 



du' 



r»L*(T,o) i T ' + i 

■ I — In ■ 

Jo t' t' - 1 M ' 2 - 2u'cos9 + 1 

where m' = w' = i r' = t== . 

max 2 T-co.sfl' Vm' 2 -2i/«m0+1 

We see that the integral depends on r and only so that 



whereC= «. 

Choosing 9 « f « 1 it will be possible to find analyt- 
ically the integrals I\ and h- Each of them is a sum of 
the two integrals: 



exp 



*-J 

r 

+ exp\ 
= /a + la 



-k 
k 



sin/3 + sind 



sin(J3 + 8) 
sin/3 + sinO 



d/3A 



sin(J3 + 8) 



cos z (J3 + 6)d/3 ■■ 



(B.4) 



dN l ° k 2 



dt & 



Thus 



dN" dN\ s 



(A.6) 



(A.7) 



So that for a given time the ratio is proportional to the 
distance R measured in units of A. It can be verified 
that this relation holds for any scattering process in a 
uniform medium, not only for isotropy. 



Appendix B. Derivation of an analytical expression 
of the angular distribution of the first 
generation dM^/dO for small angles 6 

The formula for the angular distribution of the first 
generation integrated over time, can be rewritten in the 
following form 
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k \cos8\ 
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sin(J3 + 8) 
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where f(J3 + 6) = /(or) is the angular distribution of the 
scattering angle a, normalised to 1 when integrated over 
the full solid angle. 
For the Rayleigh scattering we obtain 
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To find the above integral for 8 <K 1 we divide the inte- 
gration region into two parts: 
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for; = 1, 2andy6n = 0,/3 n = j8 2 i = Ji-8-e,/3 22 = n-0. 
Using 8 <*c 1 we have 
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exp[ - k(l + — - 8ctg/3)]d/3 = 
^ exp(-k8tg^)d/3 (B.5) 



where we have neglected all terms oc 8 2 . Note that in the 
integration region 8 ctg/3 < | «: 1 so that our approxi- 
mation 



(1 +8ctg/3T l - 1 -8 ctg/3 



(B.6) 



is justified. For the same reason the exponential index 
k8tg\ < k0 \ <sc 1 so that 



(l-k8tg^d/3 = (B.7) 
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^e- k {ii-8-e + 2k8ln-) 
The integral I\ 2 can be found easily: 
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Thus 



h = hi +I12 = 
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= e 



^n-28-2e + 2kd(2ln^ + 1) 



When calculating I 2 we introduce an angle % = n - 8 - 
/3 to take advantage of x being small in the integration 
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region of h- 
We have 



h\ = J exp\y 
Jo 



)dx 



2k6 



dx = e~ k 2k6 
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where Ce u — 0.577 is the Euler constant. 
The last integral, I22, equals 
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hi = ^ exp{-k^t(\- X 2 )d X ^ (B.ll) 
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L Jus. y Jus. y 



The second integral is of the order of e 2 , so that /22 = hi 
and 



/ 2 = e^l^etC^ - /n| + /n(ifc0) - 1] + 2e 



(B.12) 



Thus 



/1 + /2 = - 20 + 4fe0[Zn(fe0) + C Eu - i]) 

(B.13) 



and finally 

<Mff 9faT* 



rf0 



64fl 2 



(B.14) 



As should be expected, the terms depending on e have 
cancelled in the final result. 

The ratio of the scattered light arriving within a cer- 
tain angle f to the point source direction to that non- 
scattered (direct) equals 

Af* 9tt I 4kf 2 \, „ , „ , 1 1) 
-^,(*£) + C a -l-— ]J 

(B.15) 
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Appendix C. Towards an analytical solution of 
the second generation in a uniform 
medium 

Here we present results of some analytical integra- 
tions of functions F(ff, a) and G(R', a) defined in Sec- 
tion 2.2. For the Rayleigh scattering we have 



f R (a') = —(l+cos 2 a') 
16n 



(CI) 



Since cosa' = cosacosO' - sina sinff cos0' , F R (ff,a) 
(Eq. 25) can be easily found analytically 

F\ff,a)= f 2 " f\a'W= (C.2) 
Jo 



1 + \^sin 2 a + {cos 2 a - \^sin 2 a)cos 2 ff 



The calculation of a corresponding function F M (0',a) 
for the Mie scattering is a little more complicated. The 
adopted here form of f M (a') has been chosen in such 
a way as to make the integration over 0' analytically 
possible. We have (Eq. 16) 



F M (ff,a') = ff M (a')d<P' 
Jo 

= 2 f (en cos % a' + b)d<p' 
Jo 



(C3) 



The angle of the second scattering a' changes within 



limits a' min < a' < a' lnax where 



= |0' - a\ 

= ff + a if ff + a<n 

= In - (ff + a) if ff + a > n 



(C4) 



Since a' = f separates the two regions of a' where 
f M (a') has different shapes, the solution of (C.3) 
depends on the positions of a' . on the a' axis with 

< < mm 

respect to the point a' = f . There are three cases to be 
considered 



a)- aiL„ r < ?. Then, of course, a' . < \ as well and the 

/ max 2 mm 2 

function f M (a') = f^(a') applies to the whole region 
of <p' (0 < <p' < 2n). Thus we have 

F M = 2[ ^ ai(A + Bcos<f>'fdcf>' + nb 

IS r-n 
aiJ^C k I cos k <p' d<p' + nb 



(C5) 



where C k = (f)A*- k B k , A = cosacosO' and B = 
sina sinff . 

For odd k the integrals in (C.5) vanish. For even k-2m 
we have 



hr. 

Hence 



Jo 



cos Zm <f>'d<f>' = 



(2m- 1)!! 

2^n 



r 4 

F M = 2\a l Y J C 2m J 2m +nb 

m=0 



(C6) 



(C7) 
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b)- a 'min < f < a 'max- This case is a bit more com- 
plicated since there are both parts f^(a') and f 2 M (a') 
involved. The value <f>' m corresponding to a' = | equals 
cos<f>' t — ctga ■ ctgO' . 
We have 

f m = i\ f* f* t ww+ P7f(«') 

l Jo J^: 
* r r^* 

= 2 y Ck\a\ I cos k <p'd<p' + 

k=0 J° 
+ a 2 I cos k (p'd(p' 

J 01 



+ 2tt£ 



(C.8) 



Solving the integrals we finally obtain 

4 



2 c 2n {(ai - «2)[W w "k + G n (^: 

+ 2^( n") 71 "} ^ C2 " +l(fll ~ fl 2)G^(^) + for 



where 



1 ^ l2n\ 



2«\ s/n[(2« - 2/)0'J 



2n - 2i 



(C9) 



(CIO) 



/ / L sr-i ATI - 

Gn(<P*) = 7p n 2j f 



1 xi (2n + 1\ «n[(2n - 2i + 1)^] 
2n - 2i + 1 



c). a' > ? . Then, of course a' . > ? as well and for 

7 mm 2 mm 2 

all values of 0' we have that a' > |, so that f^(a') = 
f 2 M (a'). The calculation of ,F M is the same as in case a) 
but fli — > A2- Thus 

4 



F M = 2 fl 2 2 C 2m-/2m + Tib 

m=0 



(C.ll) 



The next integration over 8' (Eq. 26) defining 
G AB (R', a) can be found analytically onl>0 for G RR . The 
calculation is lengthily so that we will present here only 
the final result for the flux of the second generation 



j RR (8,t) = - — — e~- I{T',a)—r 
Jl 2\in 4 Jo V 2 



(C12) 



r 1 / 2sin 2 6' 2sin 4 6' \ 

I(T',a)= -1 _ + __ . 

Jo y\ y y 2 ' 

■{a + bcos 1 e'jsine' d& = (C.13) 

, for' 2 + l) 2 l x i 2 -x\ 

4 

z 



, *2 & 

+ flu /n 1 

x, 4r' 2 



x' 2 x\ xz 

; h fl_2'« 1- 



- 2(t' 2 + V-i^cii-i — - *' +a_i/n^j 



The coefficients a,- are the same as in (13); Xi = (t'- l) 2 , 
JC2 = (t' + l) 2 , a — 1 + ^sin 2 a, b — cos 2 a— jsin 2 a. The 
integral over x' has to be found numerically. 



where 



2 actually for G MR , G RM and G MM the result could also be found 
analytically but due to the high power of cosa (index = 8) there would 
be very many terms making such an effort worthless. 
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